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Abstract 

Currently, Eulerian flow solvers are very efficient in accurately resolving flow structures near solid boundaries. 
On the other hand, they tend to be diffusive and to dampen high-intensity vortical structures after a short 
distance away from solid boundaries. The use of high order methods and fine grids, although alleviating 
this problem, gives rise to large systems of equations that are expensive to solve. 

Lagrangian solvers, such as the regularized vortex particle method, have shown to eliminate (in practice) 
the diffusion in the wake. However, as a drawback, the modelling of solid boundaries is less accurate, more 
complex and costly than with Eulerian solvers (due to the isotropy of its computational elements). 

Given the drawbacks and advantages of both Eulerian and Lagrangian solvers the combination of both 
methods, giving rise to a hybrid solver, is advantageous. The main idea behind the hybrid solver presented 
is the following. In a region close to solid boundaries the flow is solved with an Eulerian solver, where the 
full Navier-Stokes equations are solved (possibly with an arbitrary turbulence model or DNS, the limitations 
being the computational power and the physical properties of the flow), outside of that region the flow is 
solved with a vortex particle method. 

In this work we present this hybrid scheme and verify it numerically on known 2D benchmark cases: 
dipole flow, flow around a cylinder and flow around a stalled airfoil. The success in modelling these flow 
conditions presents this hybrid approach as a promising alternative, bridging the gap between highly resolved 
and computationally intensive Eulerian CFD simulations and fast but less resolved Lagrangian simulations. 
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1. Introduction 

It is well known that the numerical simulation of advection dominated external flows is a challenging 
problem. In Eulerian flow solvers one of the biggest difficulties lies in the correct definition of boundary 
conditions, arising from the necessity to truncate the unbounded domain, see for example mmm- A 
usual expedient to circumvent this is to employ a large computational domain, placing the outflow boundaries 
sufficiently far from the region of interest, minimizing errors. Although work has been done towards the 
reduction of the size of the computational domain required for numerical stability, see for example the recent 
work in j2lj, physical accuracy still imposes strong restrictions. Another important aspect of external flows 
is that vorticity is concentrated along the wake and in thin boundary layers surrounding solid boundaries. 
Therefore, in order to make the numerical solution accurate and at the same time computationally feasible, 
adaptive schemes should be employed, e.g. [33]. In this way, high resolution is used in the support of 
vorticity whereas in the remaining regions lower resolution may be used, reducing the computational cost. 
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Nonetheless, when several moving objects are present, the continuous generation and adaptation of the 
underlying mesh is not trivial and poses several challenges, see |4411(48] and references therein. Conversely, 
the Lagrangian vortex particle method, introduced in [12] , is a mesh free method where the vorticity is carried 
by particles which are then transported by the flow field. In this way, the compact vortical structures that are 
characteristic of the flow in wakes are naturally captured by this method, which is intrinsically adaptive since 
vortex particles are placed only where vorticity exists, }43j . Additionally, boundary conditions at infinity 
are automatically satisfied. Although many developments have made this method a robust alternative 
to Eulerian flow solvers, some challenges still exist. For example, the enforcement of viscous boundary 
conditions, introduced in [55] and improved in I54j . is still a challenge on complex geometries, see mmm- 

An alternative approach, followed in this work, is to divide the computational domain into different 
sub-domains and employ solvers specifically tailored to the flow characteristics present in each region. This 
gives rise to what is commonly referred to as a hybrid solver. Within this framework it is possible to use 
Eulerian flow solvers in the vicinity of solid boundaries and a Lagrangian vortex particle solver for the wake. 
In this way, the Eulerian flow solvers are used to efficiently and accurately resolve the regions where viscous 
effects dominate the convective ones. For instance, the flow in boundary layers, which is predominantly 
anisotropic and unidirectional, is particularly suited for an Eulerian formulation where elongated cells may 
be used. On the other hand, the vortex particle method accurately and efficiently captures the wake flow. 
Moreover, by using a multi-domain formulation other advantages arise. First, since meshes exist only in the 
vicinity of solid boundaries, the movement of each body can be treated independently of all the others, thus 
greatly simplifying the implementation of moving bodies. Secondly, despite less favourable results in the 
past, [55, 55], recent results start to establish that fast particle solvers such as the Fast Multipole Method 
(FMM), [U[25], can scale to massively parallel computations more efficiently than matrix and Fast Fourier 
Transform (FFT) solvers, [3, U, 751. By employing a hybrid formulation, large scale simulations involving 
several objects, may replace a single large Eulerian solver by a set of smaller ones coupled to an efficient 
mesh free vortex particle solver where velocities are computed by a FMM. Multidomain Eulerian solvers 
are well established, nevertheless, the authors consider that the prior division of the computational domain 
into disconnected sub-domains potentially leads to a more efficient solver. These potential improvements 
in computational efficiency have several applications, namely on the simulation of large wind farms where 
large numbers of objects move with respect to each other, different flow scales interact and the accurate 
propagation of the wake is fundamental to determine the optimal placement of wind turbines. 

1.1. Literature review of hybrid grid-particle flow solvers 

Although hybrid grid-particle methods are not a popular approach for solving flow problems, several 
formulations have been proposed in the past. 

Particle-grid methods have been formally introduced for the first time by Cottet, [15] . where a (it, ui) 
formulation of the Navier-Stokes equations is used in the whole domain. In that work, the Eulerian and 
Lagrangian domains partially overlap and the interface conditions are obtained by an alternating Schwarz 
method. This approach was extended by Ould-Salihi et al., ED, to a (it, p) formulation of the Navier-Stokes 
equations. A more detailed description of this hybrid formulation is presented by Cottet and Koumoutsakos, 
m > both for the (it, ui) and a (u, p ) forms of the Navier-Stokes equations. 

One of the few formulations without overlap is the one by Guermond et al., 1221EQ1. In this approach the 
(ip, to) form of the Navier-Stokes equations is solved in the Eulerian domain by a finite differences method. 
The transmission condition is obtained by assuming that viscosity is small in comparison to the convective 
term. 

Huberson and Voutsinas present an interesting overview of particle-grid methods in their overview article, 

CHI- 

More recently, two approaches related to the work presented here have been proposed that do not rely 
on the classical alternating Schwartz procedure. Daeninck, m, presents a hybrid formulation where the 
particles cover the whole fluid domain and completely overlap the Eulerian one. The Eulerian solution is 
then used as a correction to the particle solver in the vicinity of the solid boundaries. Stock, Gharakhani 
and Stone, [66] . proposed an improved version of this method in order to avoid interpolation errors within 
the viscous boundary layer, which for high Reynolds numbers can be large. 
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Oxley, [52], and more recently Papadakis and Voutsinas, [55], present a 2D hybrid compressible Euler 
method for transonic rotorcraft applications using a complete overlap and an iteration procedure to eliminate 
projection errors between the particles and the grid and ensure compatibility between the solutions in the 
two computational domains. 

With a focus on computer graphics animation, Golas et ah, (25], present a partially overlapping domain 
decomposition approach based on a (u, to) formulation of the Navier-Stokes equations. 

1.2. Outline 

As mentioned before, the goal of this work is to present a hybrid grid-particle flow solver. The main 
principle behind this hybrid approach is to decompose the fluid domain into a set of sub-domains where the 
most suitable flow solver is applied. In the vicinity of the boundary regions an Eulerian solver will be used, 
whereas the wake will be solved by Lagrangian particles. 

Therefore, the outline of this paper is as follows. In Section [2] we present the hybrid grid-particle solver, 
starting with a brief introduction to its key ingredients in Section [2~T] This is followed by a discussion of 
the Lagrangian and Eulerian solvers in Section |2.2| and Section |2.3[ respectively. In Section |2.4[ we present 
in detail the coupling strategy between the two flow solvers. After introducing our approach, we apply it to 
different test cases in Section [3] Finally, in Section [4] the conclusions and further outlook of this work are 
discussed. 

2. The hybrid flow solver 

2.1. Hybrid solver introduction 

As mentioned before, the objective of a hybrid solver is to use the most suitable numerical method in each 
region. In the case of the grid-particle approach discussed in this article, Eulerian grid solvers are used to 
discretize the flow equations in the vicinity of solid boundaries and a Lagrangian vortex particle formulation 
approximates the flow in the wake region. Most of the hybrid approaches discussed in Section [O] are based 
on domain decomposition with partial overlap, as depicted in Figure [T] Another common characteristic of 
most of these formulations is the use of Schwarz’s alternating method to ensure compatibility between the 
two solvers. 



Figure 1: Standard geometrical decomposition of the flow domain with partial overlap. The Eulerian domain, 1 1 . is a 
neighbourhood of the solid boundaries and the Lagrangian domain, Q L , is a neighbourhood of infinity. Both domains overlap, 
such that fis n 7^ 0. 

Although the Schwarz alternating method is an excellent technique to obtain a very accurate match 
between the solutions of the two solvers in the overlap region, D flj,, it requires iterations. The need to 
iterate makes this approach computationally expensive. For this reason, in this work we extend the ideas 
pioneered by Daeninck, m, Stock, Gharakhani and Stone, [fifil . and propose a vorticity conserving hybrid 
grid-particle solver that does not require an iterative process to ensure compatibility. 

In this work, the Lagrangian domain, completely overlaps the Eulerian one, as can be seen in 
Figure [2jr. Instead of focusing on exactly matching the solution of both solvers in the overlap region, this 
approach follows another route. The Lagrangian vortex particle method is used to obtain the flow solution 
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Figure 2: Representation of the different domains present in the domain decomposition strategy: presented in this work (a) 
Hybrid domain, (b) Eulerian domain and (c) Lagrangian domain. The Lagrangian domain solves the fluid flow in the whole 
domain and the Eulerian one plays the role of a near-wall correction. 


in the whole domain, Figure |2j:, without resolving the boundary layers near the solid walls. The Eulerian 
solver is then used to correct it in the near-wall region, Figure [2 Jd. 

To evolve from time instant t n to t n + 1 , the coupling procedure of the hybrid grid-particle solver discussed 
in this work can be summarised in the following steps: 

1. Evolve Lagrangian solution: Evolve the vortex particles from time instant t n to t n+ \ neglecting 
vortex generation at the solid boundaries, following a standard procedure as presented, for example, 
by Cottet and Koumoutsakos in m- 

2. Determine Eulerian boundary conditions: Use the Lagrangian solution at time instant t n .pi to 
compute the Dirichlet boundary conditions for the velocity field in the Eulerian domain. 

3. Evolve Eulerian solution: Evolve the Eulerian solution from time instant t n to t n+ 1 (possibly using 
kE sub steps in order to ensure stability of the solution), using any standard Navier-Stokes grid solver 
based on the (u,p) formulation. 

4. Correct Lagrangian solution: Use the Eulerian solution at time instant t rl+ \ to correct the La¬ 
grangian one in the near wall region, n Ql . 

These steps are performed in a loop, as depicted in the flowchart presented in Figure [3] 


start 



Figure 3: Flowchart of the coupling procedure to evolve the hybrid grid-particle solver from time instant t n to t n + l. 

In order to introduce the hybrid grid-particle solver we will first discuss the Lagrangian vortex particle 
method used in this work, followed by the presentation of the finite element discretization employed in the 
Eulerian domain. 
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2.2. Lagrangian particle solver 

Since this work is intended for a wider audience, potentially unfamiliar with the vortex particle method 
(VPM), in this section we give a brief introduction to this topic, emphasising the specific details of the 
implementation used in this work. For a more detailed overview of vortex methods we recommend the 
extensive reference book by Cottet and Koumoutsakos, PH; and the briefer reviews by Raviart, mi and 
Winckelmans, m- Since this work focuses on two dimensional flows, we will restrict our analysis to two 
dimensional domains. 

On a simply connected two dimensional unbounded domain, Q = R 2 , and for time interval t £ (0, T], 
the system of the Navier-Stokes (NS) equations describing incompressible fluid dynamics in the (u,u>) 
formulation in the absence of external forces is given by: 


duj 

— + (u ■ V) to — v\I 2 uj = 0 , 

in fl , 


(1) 

V • u = 0 , 

in Q , 


(2) 

V x u = uj . 

in Q , 


(3) 

u(x,t) = UJ 0 (x) , 

in fl and for 

t = 0, 

(4) 

lim u(x,t) = U 00 . 

in Q and for 

ie]0,T], 

(5) 

lim oj(x, t) = 0 . 

in fl and for 

f €]0,T], 

(6) 


If we consider inviscid flow, v = 0, (JT]) reduces to: 


Doj 

~Dt 


dui 

~dt 


+ (u • V) uj 


0 , 


(7) 


where ^ is the material derivative. This equation simply states that vorticity is transported along the 
velocity field lines, see Figure [4] This points to a Lagrangian formulation where fluid parcels containing 
vorticity are traced along field lines. 



Figure 4: Fluid parcels containing vorticity are traced along field lines. Where the positions of fluid parcels at time instant 
t = to are represented by [•] and the positions of fluid parcels at a time instant t > to are represented by [•]. 

The vortex particle method stems from this Lagrangian formulation. For incompressible and inviscid 
flows, the vortex particle method assumes that fluid parcels convect without deformation, therefore the 
position, x p , of each particle, p , evolves according to the following ordinary differential equation: 

= u(x p ,t), (8) 

where u(x, t) is the velocity field associated with the vorticity distribution uj(x, t) due to the vortex particles. 
Given the incompressibility constraint, ([2|, and the definition of vorticity, ([3]), it is possible to show that 
the velocity, w, and vorticity, ui, are related by the following Poisson’s equation: 

V 2 « = V x (ue z ) , (9) 


with boundary conditions as given in <©• 
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Green’s function formulation can be used to directly compute the velocity field from the vorticity field 
by noting that the velocity, u$, associated with a Dirac-5 distribution of vorticity, u>s, located at x p satisfies 

V~u s = V x (oj s e z ) := V x (S(x - x p )e z ) . 

For homogeneous boundary conditions, us is given by the expression: 

1 x — x p 

u 5 ~ x 77 rry e z ■ 

2tt ||£c - x p \\ 2 

Due to the linearity of Poisson’s equation, a given vorticity distribution, w( x), together with the boundary 
conditions <[5j) generates the following velocity field: 


uix) = -hJ a ^Xp xeMi)dS:+u ^- (10) 

If we associate to each vortex particle, p , a vorticity distribution, u) p (x), the total vorticity distribution, 
ui(x), will be given by 


u(x) =' 52 u} p (x) , 


( 11 ) 


and the total velocity field will be, replacing (11) into (10) 


1 


t(x) = -V 

V ’ 2tt ^ 


X - x 


x e z ujp(x )di + Uao . 


( 12 ) 


In the vortex particle method the vorticity distribution associated with each particle, u> p (x), is given by 

w p (®) :=C<7(||*-*p||)r p , (13) 

where T p := u(x)dx is the circulation contained in the vortical fluid element associated with each 

vortex particle, and C<r( r ) := :7yy2C(y) is referred to in the literature as regularised kernel, see for example 
mm- The parameter a is the core size associated to the size of each vortex particle. Typically T p is 
computed by numerical integration using some kind of quadrature rule. The most common one, and used 
in this work, is to use lowest order Gauss-Lobatto quadrature, 


T p = (jj(x p )h 2 . 


Under these conditions, if we substitute (|T3 1 into ( |12| ) we get 

9a(\\X ~ 




\x — X Tl 


(x Xp ) x & z r p , 


(14) 


(15) 


where g a (r) := g(^). Several options exist for the pair of functions (gX), see for example [T7, 72, for the 
more standard choices, or mm for a more extensive discussion of Gaussian kernels and (i(H [65], El] for 
arbitrary order algebraic kernels. In this work we use the more standard Gaussian kernels of order two: 


g(p) := 1 — < 


and 


C(p) ■■= 


(16) 


For the evaluation of the velocity held by (15) we use a Fast Multipole Method for a fast and grid-free 
computation, see Goude and Engblom, [25], and Engblom, [22], for further details. 

The convergence proof for vortex particle methods, first introduced by Hald in [31] . extended simultane¬ 
ously by Beale and Majda in [3] and Cottet in m , is presented in full detail by Cottet and Koumoutsakos in 
m- This proof establishes the mathematical justification for the success of this method. On the other hand, 


6 










it also establishes an important upper bound on the inter-particle distance, h, required for the convergence 
of the method, namely 

h<Ca 1+s with s, C > 0 . (17) 

Due to the nature of fluid flow, it is known that there exists a time T after which © is no longer satisfied 
and the method no longer converges. Several techniques have been introduced to preclude the inter-particle 
distance from increasing above the bound ( |17| ). Iterative circulation processing, [S], rezoning , M, and 
remeshing , |36| , have been three of the most popular approaches. In this work we opted for remeshing for 
its computational efficiency. This method introduces an underlying structured mesh with spacing equal to 
a target inter-particle distance, h, satisfying the overlap criterion ©■ At regular time steps the vortex 
particles are reinitialised into this grid. This reinitialisation consists of generating a new set of particles on 
the regular mesh by means of conservative interpolatory kernels with compact support, W, according to the 
expression 

I7 W = E w « ew - < d >vr-y°v A ) r ? d - ( 18 ) 

P&Qq 

where the particle, p, located at x = [x p ,y p ] has an associated circulation T p and Q q is the set of indices 
associated to the particles that lie in the support of W (x q ew — x, y q ew — y). Typically, this higher dimensional 
interpolatory kernels are constructed as the tensor product of one dimensional ones, 

W&rj) :=W(OW( 77 ). 


Many options exist for the definition of W. each with different orders of accuracy and properties, see for 
example mm for a detailed discussion. In this work we use the popular M' A interpolating kernel introduced 
by Monaghan in [48| and given by the following expression 


K(p) ~ 


l-fp 2 + !|p| 3 if \p\ < 1 

\{2~\p\f{l-\p\) if 1 < |p| < 2 
0 if |p| >2, 


(19) 


and IV(£) = M' 4 (j). This kernel is known for its good combination of smoothness and accuracy, as pointed 
out by Cottet et al. in [TS] . 

For viscous flows, v > 0, viscosity is not only transported along the velocity field lines but also diffused. 
Therefore ([8]) is extended into 


i da? p 
d t 
dw 
d t ~ 


- u(x p , t), 
^V 2 w. 


( 20 ) 

( 21 ) 


This characteristic of fluid flow led to the introduction of the method of viscous splitting by Chorin in 1121 . 
Instead of solving the coupled system of convection (20) and diffusion (21), this method decouples the two 
equations. This is done by introducing two sub-steps. In the first one, particles move along the field lines, 
according to (20), as discussed before. In the second one, the diffusion effects are computed according to a 
discrete version of (21). The most popular approaches for the discretization of the diffusion term, (21), are 
the particle strength exchange (PSE) method, (20] . the vortex redistribution method (VRM), [63] . and the 
random vortex method (RVM), [T21 . For a detailed discussion of the different particle diffusion schemes we 
recommend the shorter overview by Barba, [2], or the more extensive works by Barba, jT], and Cottet and 
Koumoutsakos, m ■ Although the VRM has the advantage of maintaining the Lagrangian character of the 
vortex particle method, since it does not require an underlying mesh, it comes at the cost of solving at each 
time step an underdetermined system of equations per particle. In order to improve the efficiency of this 
method, Wee and Ghoniem, [65], proposed to redistribute the particle strengths into target particles placed 
on a regular mesh. Although very efficient, this method imposes restrictive bounds on the size of the time 
steps, particularly on their minimum size. To circumvent this limitation, we opted to employ the method 
of Wee and Ghoniem within the context of the Particle-Mesh approach used by Sbalzarini et al. [62]. This 
formulation consists of two steps: 
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1. Remesh the particles into a regular mesh, using (18) together with the interpolation kernel, (21). 

2. Diffuse the vorticity by redistributing the particle’s circulation according to 




1,3 


1,3 


1 - 


AvM Av 2 At 2 


h 2 


h 4 


+ (Tl hj + rf +1)i + !*•_, + r^. +1 ) -ry - 


( vA t 2v 2 At 2 

h 4 


+ ( r f-iy-i + r ?-i j+i + r 




V h 2 

~\k \ 

■ i+l,j+l) 


2 At 2 


h 4 


( 22 ) 


In this expression, as before, h is the inter-particle distance in the regular mesh, the indices i,j correspond 
to the particles’ locations on the mesh, Xij = [ih,jh], At is the time step size and the index k corresponds 
to the time instant t & = kAt. This particle circulation update corresponds to the diffusive redistribution 
scheme presented in jB5j for the special case of regularly distributed particles. It is important to note that 
this update also satisfies the redistribution formulae given in [55] and is equivalent to a PSE method, DU- 
Moreover, (22), corresponds to a finite difference in time solution of the diffusion equation on a mesh, similar 


to what is done in [62] . 

For viscous flows in the presence of solid boundaries, two additional issues need to be taken into account: 
(i) the enforcement of the no-slip boundary conditions and (ii) the vorticity flux at the solid boundaries. 
Currently, the standard approach is based on the vorticity boundary conditions pioneered by Lighthill, [45] , 
and then extended to the vortex particle method by Cottet, BE and Koumoutsakos et al. in [IWFuT] . For 
a more detailed discussion of this and other approaches, the reader is directed to the book by Cottet and 
Koumoutsakos, D2I- On what follows, we give a brief description of this method and its differences in the 
context of the hybrid flow solver. 

In summary, this formulation consists of combining the vortex particle method with the boundary element 
method (vortex sheet panels). Initially, at the beginning of time instant t n = nAt, it is assumed that a 
vorticity field satisfying the no-slip boundary condition exists. The first sub step of this algorithm consists 
in computing the advection-diffusion evolution of the vortex particles, (20) and (21), as explained above, 


therefore without enforcing the boundary conditions. By the end of this first sub step, since no boundary 
conditions are enforced, a spurious slip velocity exists. In the second sub step the vortex sheet strength that 
cancels this slip velocity is computed. This in turn is related to the vorticity flux at the solid boundary. 
This method is represented graphically in Figure [5] left. 


start start 



Figure 5: Flowcharts of the viscous Lagrangian vortex particle method (left) and hybrid method (right) in the presence of solid 
boundaries. 


In the hybrid method presented in this work, the computation of the vorticity flux is intrinsically taken 
into account in the Eulerian sub-domain and transferred to the Lagrangian one in the correction sub-step 
(see Figure [dj) . Therefore, in this work, the computation of the vorticity flux is replaced by the correction 
step and the computation of the vortex sheet strength, 7 (s), is performed at the end in order to enforce the 
no-slip velocity boundary conditions, as depicted in Figure [5j right. 
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In this work we focus only on enforcing the no-slip boundary conditions since, due to the linked boundary 
conditions , there is a direct relation between the no-slip and the no-through boundary conditions. This 
relation is detailed by Koumoutsakos, [35], and Shiels, [64 1. 

The vortex sheet strength is obtained from the following Fredholm equation of the second kind: 

“ x ( s ')l)] TOOds' = Uelip ' s , (23) 

where T, w is the solid boundary, u s u p is the slip velocity and s is a unit vector tangent to the solid boundary 
in the direction of integration. These equations were discretized using the boundary element method (panel 
method), resulting in: 

^ 2 “^ ~ 9 [ lo g(|*(s) - ®(s')|)] 7fc(s') ds ' = u sli P ■ s for fc = l,...,JV, (24) 

n k= 1 ^ ak n 


where 7 ^ is the vortex sheet distribution associated to the panel k and <ik is the boundary segment associated 


to the panel k. The drawback of this formulation is that the system of equations (24) is singular, and therefore 


allows an infinite number of solutions. An additional equation that limits the solution space to a unique 
solution is one prescribing the total strength of the vortex sheet: 


N r 

^2 7fe(s)ds = r 7 . 
fc =1 


(25) 


The total strength, T 7 , can be obtained from Kelvin’s circulation theorem which states that the total 
circulation must be conserved in time. As will be seen in more detail in Section [274] the total circulation of 
the vortex sheet can be determined from the coupling between the Eulerian and the Lagrangian sub-domains. 

By adding this additional equation, (25), the resulting system of equations will be overdetermined since 
there will be N unknowns but N + 1 equations. Several options exist to solve this problem, see for example 
Cottet and Koumoutsakos, In this work we solve this overdetermined system of equations using a 

Least-Squares approach. 

It is worth to mention a more robust formulation proposed by Koumoutsakos and Leonard, EH, consisting 
of eliminating the first term in the spectral decomposition of the kernel of the system of equations. By doing 
this, the singularity is removed and the system of equations is well-posed, in the future we intend to explore 
this route. 

By following the steps outlined above, at each time step the velocity held in the Lagrangian domain will 
be given by: 

U = U OQ T T Ury , (26) 

where U oo is the free stream velocity, u u is the velocity associated to the vortex particle distribution and 
u 7 is the velocity induced by the vortex sheet at the solid boundary. 


2.3. Eulerian solver 

In this section we introduce the Eulerian how solver used in the hybrid method presented in this work. 
The most common approaches to discretize the Navier-Stokes equations are the Finite Volume Method, the 
Finite Differences Method and the Finite Element method. We have opted for a Finite Element discretization 
for its geometrical flexibility and fop-refinement capabilities, which we intend to use in future developments 
of this work. 

Here, as opposed to Section 2.2 we consider the (u,p) formulation of the Navier-Stokes equations in a 
two dimensional bounded domain, f2 C K 2 , in the absence of external forces, 


— + (m-V)m-V-ct = 0, 

in Q , 

(27) 

V ■ u = 0, 

in Q , 

(28) 

u(x,t) = Uq(x) , 

in fl and for t = 0 , 

(29) 

u(x, t) = Ub(x, t) . 

in Oil and for t G] 0 , T ], 

(30) 
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where u denotes the velocity, u q corresponds to the initial velocity condition, u & is the velocity at the 
boundary and a is the Cauchy stress tensor defined as 

<j(u, p ) := 2 i/e(u) — p I, 

where v is the kinematic viscosity, p is the pressure, I is the identity matrix and e(u) := |(Vtt + Vid). 

To solve the Navier-Stokes equations with the Finite Element Method (FEM) it is necessary to construct 
the associated weak form and to select the appropriate function spaces for the velocity, u , and pressure, p, 
in order to ensure stability. 

It is well known that suitable function spaces must satisfy the Ladyzhenskaya-Babuska-Brezzi (LBB) 
compatibility condition, see Brezzi and Fortin, [5]. Such a pair of function spaces is the Taylor-Hood family, 
m- This family of function spaces consists in C° Lagrange interpolants of polynomial order two (V/,) for 
the velocity and polynomial order one ( Qh ) for pressure. 

For the construction of the weak form we use the Incremental Pressure Correction Scheme (IPCS). This 
method, introduced by Goda, [23], is an improvement to Chorin’s projection scheme, and consists of 


computing a tentative velocity, u *, by neglecting the pressure in (27) and then correcting it by determining 


the pressure field that gives rise to a divergence free velocity field. To advance the solution in time we use 
a forward Euler scheme. This algorithm can be summarised in the following three steps: 

1. Compute the tentative velocity: At time instant t n = nAt , find the approximate tentative velocity 
u* h G Vh such that: 


At 


-,v) + « • V<-\ v) + {o{ul-\pl~'),e{v)) 


(31) 


+ (p% L fi, v) d n - (vh • (Vu ft 2 )\ v) d n = 0 Vi> e V h , 


2 . 


where h is the unit vector normal to the boundary, dfl, and u h 2 := Uh —• The Dirichlet boundary 
conditions for Uh, (30), are also applied to u £ in this step. In the hybrid method, these correspond to 
no-slip conditions at the solid boundaries and at the external boundary of the Eulerian sub-domain, 
Ed in Figure [2] the velocity is obtained from the Lagrangian sub-domain. 

Compute the pressure: The pressure is obtained by finding p ^ £ Qh such that: 


(VpJJ.Vg) = <Vpr\Vg> - ^(V -u* h ,q) 


V? £ Qh ■ 


(32) 


Where weak homogeneous Neumann boundary conditions are used for the pressure, see Guermond et 
al. [2S] for a discussion on pressure boundary conditions in the context of projection schemes. 

3. Determine the corrected velocity: The corrected velocity field is obtained by finding u r t \ £ Vh 
such that: 

( u h, v ) = { u h, v ) - At (^iPh -Ph -1 )^) v£V h . (33) 

The Dirichlet boundary conditions for velocity are taken into account at in step 1. 


This solver was implemented using the finite element library FEniCS, see Logg and Wells, 
details. 


for more 


2-4- Hybrid solver 


As seen in the introduction to the hybrid flow solver, Section 2.1 the main idea behind the hybrid 


solver discussed in this work is the following. Divide the computational domain into two sub-domains: one 
Eulerian and one Lagrangian, see Figure [2] The Lagrangian solver covers the whole computational domain 
and the Eulerian one exists only in the vicinity of solid walls, being capable of capturing near-wall effects and 
serves as a correction for the vortex particles in that region. This, in turn, results in an evolution algorithm 
comprised of four steps: evolve Lagrangian, determine Eulerian boundary conditions, evolve Eulerian and 
correct Lagrangian, see Figure [3] Therefore, in this section, we will discuss each step of the evolution 
algorithm of the hybrid method. 
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2-4-1- Evolve Lagrangian 

The first step of the algorithm to advance the hybrid solver solution from time instant t n to t n +i = 
t n + A tL is the evolution of the Lagrangian sub-domain. We assume that at the start of the algorithm the 
velocity field in the overlap region is nearly identical in the two sub-domains. Therefore, at the start of this 
step the velocity field in the Lagrangian sub-domain, (26), satisfies the no-slip boundary conditions. After 


evolving the Lagrangian sub-domain the vortex particles are remeshed into a regular grid, using (18), as 
discussed before. 


Under these conditions, the Lagrangian sub-domain is evolved as described in Section 2.2 Once more 


we note that, in the context of the hybrid solver, the evolution of the Lagrangian solver does not explicitly 
include the generation of vorticity at the solid boundary. For this reason, at the end of this step, the 
Lagrangian solution will not fully resolve the near-wall region. Nevertheless, as discussed by Daeninck, m, 
and Stock, Gharakhani and Stone, jaB], this error will not affect the accuracy of the induced velocity field 
in the outer boundary of the Eulerian sub-domain, Ed, Figure [7] right. 

2-4-2. Determine Eulerian boundary conditions 


Once the Lagrangian solution is evolved, it is possible to use (26) to compute the velocity field, u b , at 
the outer boundary of the Eulerian sub-domain, Ed, at time instant t- n +i =t n + A tL- If the velocity field at 
the outer boundary is required at any other time instant t £ [t n ,t n +i] a linear interpolation in time is used, 


u b (t) = u b (t n ) + 


t ty j 

Air 


\u b (t n -\-l) n,b{t-n)\ t £ . 


(34) 


Note that higher order interpolation could be used. 


2-4-3. Evolve Eulerian 

Due to the different nature of the solvers in the two sub-domains, different time step constraints are 
required. Typically, the time step required in the Lagrangian sub-domain, A t^, will be larger than the one in 
the Eulerian sub-domain, A tE < A t^. For this reason, and in order to improve the computational efficiency, 
each solver is allowed to advance in time according to different time steps, subject to: A = kE^ts with 
kE £ N, see Figure [6] 


A tE 



Figure 6: Eulerian multi-stepping to match the Lagrangian time step. The figure shows A tL = kE&tE with (kE =4) Eulerian 
sub-steps to time march from t n to t n +i- 


Therefore, to advance the Eulerian sub-domain from time instant t n to t n + 1 , kE Eulerian evolution 
sub-steps associated to the time instants 


fk — j. 


kAt-i 


k = 0,...,k E 


are performed. The velocity boundary conditions at the intermediate time instants are given by (34). 
These Eulerian evolution sub-steps are performed according to what was described in Section [2~TT| 
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2.4-4- Correct Lagrangian 

The final step of the algorithm for evolution of the hybrid solver is the correction step. As mentioned 
before, in this step the vorticity field computed in the Eulerian sub-domain is transferred to the Lagrangian 
one. This step replaces the vorticity flux step in pure Lagrangian vortex particle solvers. This step is 
fundamental for an accurate coupling between the two sub-domains, therefore it is crucial that this transfer 
of vorticity satisfies conservation of circulation. 

In this work, as in the work of Stock, Gharakhani and Stone, j66j . the interpolation (or correction) region, 
tlf. is a subset of the Eulerian sub-domain, f tj C He, see Figure[TJ left. The interpolation region is defined 
such that its outer boundary, E 0 , is at a distance d^dry of the outer boundary of the Eulerian domain, E^, 
and its inner boundary, E i; is at a distance d sur f of the solid wall, E^,, see Figure [7J right. This is done 
firstly to ensure that the high gradients of vorticity near the solid boundary do not introduce interpolation 
errors in the Lagrangian domain and secondly to minimise the errors in the Lagrangian sub-domain due to 
small discrepancies in the boundary conditions of the Eulerian sub-domain. 



„ o O 

-o-FTtV 

I„ O „ _ CW -£ 7^ 


o o 

Qj 

tX 





body 


Figure 7: Definition of the Eulerian, Q, and Lagrangian, sub-domains and interpolation domain flj where the Lagrangian 
solution is corrected, with boundary dfti = E, U E 0 . Left: Representation of the whole domain. Right: Detail of the sub- 
domains close to the boundary. 


The algorithm used to correct the Lagrangian solution in the vicinity of the solid walls is summarised 
as: generate new particles, assign circulations and ensure no-slip conditions. 

The first step in the correction of the Lagrangian solution is the generation of new particles in the 
interpolation region, Qj. We start by identifying the particles that lie inside the region bounded by the 
outer boundary of the interpolation region, E 0 , red region in Figure [8] left. These particles are then removed, 
Figure [8] centre. We do not use the region Qi as a selection region since, during the evolution of the 
Lagrangian solution, some particles may move to the region bounded by E* and E„ and therefore would not 
be removed. Finally, the interpolation region is filled with zero circulation particles regularly distributed, 
Figure [8] right. 




Figure 8: Correction of Lagrangian solution in the vicinity of a solid object k : generation of new particles. Left: Particles in 
the region enclosed by the boundary Ejj are selected, centre: These particles are removed. Right: The interpolation region, 
£7/, is covered with new particles with zero circulation and regularly distributed. 


The second step corresponds to the assignment of circulation to the newly generated vortex particles. 
This follows the same procedure as outlined in Section 2.2, namely (14). Now, the vorticity field is the one 
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obtained from the Eulerian finite element solver, u>h(x, t) = V x Uh(x 7 1). In general the Eulerian sub-domain 
will be discretized with an unstructured grid. Evaluating the discrete vorticity field at the location of each 
new vortex particle is computationally expensive when unstructured grids are used. For this reason, we use 
an intermediate fixed structured grid. The vorticity field is initially interpolated into the structured grid 
and then from the structured grid into the vortex particles. 

The structured grid is constructed such that is covers the unstructured one, see Figure [To| left. Since the 
two grids are fixed to each other, the weights, Wij. that define the matrix W representing the interpolation of 
vorticity from the unstructured grid, u>h, onto the structured one, Cjh, are constant throughout the simulation 
and therefore need to be computed only once, reducing considerably the computational cost. We can then 
compute the vorticity in the structured grid, see Figure [9j by using 

&h,i = ^ ^ WijW h ,j • 


w 



Figure 9: Interpolation of the vorticity lj from the unstructured Eulerian grid into a fixed structured grid. 


The computation of the vorticity field at the particles’ locations is done by first transforming their coor- 
dinates into the reference frame of the structured grid, see Figure[lO]middle. This allows a fast identification 
of the cell containing each vortex particle. The computation of the vorticity at each particle is obtained by 
linear interpolation, 

4 

{Xp ) ^ ^ ^h,q i 

9=1 


where q € {1,2,3,4} denotes the four nodes of the cell in which 


Figure 10 right, (bh 7 q represents the associated values of vorticity in 


of the bilinear interpolation given as usual by, 


the particle p is located, as shown in 
those vertices and W pq are the weights 


W p i 
W p2 

Wp 3 

W p4 


(d' x , p - h)(d' VtP - h) 

h 2 

dx, P {dyp h) 

h 2 ’ 


h 2 ’ 

dy 7 p(d x ,p /z) 

h~ 2 


where d' = [d' x p ,d' yp ] is the coordinates of the particle p in the reference frame of the cell, see Figure 10 
right, and h is the cell spacing of the structured grid in both x and y directions. We use the same grid spacing 
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in the structured grid as in the remeshing of the vortex particles. Finally, the circulations are assigned to 


each particle, using (14), see Figure 10 



Figure 10: Algorithm to assign circulation strengths to the particles from the vorticity distribution in the structured grid. Left: 
Vortex particles overlapped to the structured grid and Eulerian unstructured grid, centre: Vortex particle, p , in the reference 
frame of the structured grid is represented by d p . Right: Vortex particles in the reference frame of a cell of the structured grid 
are represented by d' p . 


uh 2 



Figure 11: Assignment of circulation into the new blobs in the Lagrangian domain. 


The third and final step consists in ensuring the no-slip boundary condition since, at this stage, there 
exists a slip velocity at the solid walls. As discussed in Section |2.2[ we need to compute the new vortex 
sheet strength that satisfies the no-slip boundary conditions. Because we use a hybrid approach, we do not 
need to compute the vorticity flux at the boundary, since the assignment of circulations to the new particles 
automatically takes this into account. It was also seen in Section 2.2 that the new vortex sheet can be 


obtained by solving a discrete version of the Fredholm integral equation of the second kind, (24) and the 


singularity of this equation circumvented by prescribing the total circulation of the vortex sheet, r 7 , as an 
additional equation and solving the overdetermined system in a Least-Squares sense. For bodies rotating 
at a constant speed, Kelvin’s circulation theorem states that the total circulation of our system remains 
constant in time. Therefore we have: 


r 7 (f) +r L (t) = constant, (35) 

where Tl := ^ p Fp, is the total circulation contained in the vortex particles. For the case of bodies with a 
time varying rotation a similar expression can be deduced, see S3- Since the Eulerian and the Lagrangian 
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solvers must coincide in the correction region, f1/, we have that: 


V L(t)\n I = u h (x,t)dx, 
J Qt 


( 36 ) 


where Fx, (t) denotes the expected total circulation contained in the vortex particles that lie in the 

correction region. Since the step where the circulation is assigned to the particles in the correction region 


relies on an approximate quadrature formula, (14), there will exist an error, <51^, therefore: 


T L (t)\ n = +6T L (t)= uj h {x,t)dx , 

Jih 

where r ^(t) is the actual total circulation of the vortex particles that lie in fl/. In order to exactly 

Oj 

conserve circulation, (36), the circulation of the newly generated vortex particles is corrected using: 

sr L 


-'corrected 


— r„ + 


N ni ’ 


where Nq 2 is the number of particles that lie inside the correction region, Qj. 
To finally compute the total circulation of the vortex sheet notice that: 


T 1 (t)+T L (t)\ QE = f u h (x,t)dx. 

J Qe 

Therefore the total circulation of the vortex sheet can be simply obtained by the expression: 

r ■y(t)= f u h (x,t)dx-T L (t)\ nE . 

J Qe 

It is important to note that if the Lagrangian solver could fully capture the boundary layer, then the total 
circulation of the vortex sheet, T 7 , would be zero, although each vortex panel would not be zero. We see, 
then, that the vortex sheet plays not only the fundamental role of ensuring the no-slip boundary conditions, 
but also of representing the boundary layer. By following this approach, the hybrid solver exactly conserves 
circulation, greatly improving the results, see m for further details. 


3. Numerical benchmark cases 

On what follows we present the results for the solution of different flow problem using the hybrid solver 
presented in this work. In order to adequately compare the full FEM solution to the hybrid solution we 
have constructed the meshes such that they are identical in the Eulerian sub-domain. 


3.1. Dipole in unbounded flow 


The first test case is the evolution of a vortex dipole in an unbounded domain. With this test case we 
intend to investigate how the flow solution produced by the hybrid solver is perturbed as it traverses the 
Eulerian sub-domain. In order to do this, we used as initial condition the Clercx-Bruneau dipole, [13], with 
a positive monopole at (xi,yi) = (—1.0, 0.1) and negative monopole at (cc 2 , J/ 2 ) = (—1.0,—0.1), each having 
a core radius R = 0.1 and characteristic vorticity magnitude ui e = 299.528385375226 as given by Renac et. 


ah, [58], 


w(x,y, 0) = w e 




(37) 


where rf = (x — Xi) 2 + (y — yfl 2 . The Eulerian sub-domain is defined as He = [—0.25, 0.25] x [—0.5, 0.5], 
meaning that the dipole is initially placed to its left, as depicted in Figure [l2| 
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Figure 12: The domain decomposition for the Clercx-Bruneau dipole convection problem, with the positive pole located at 
p_l_ = (# 1 , 3 / 1 ) = (— 1 , 0 . 1 ) and negative pole located at p— = (# 2 , 3 / 2 ) = (—1 , — 0 . 1 ) . ( Not to scale) 


We ran a simulation with the hybrid solver, using the parameters presented in Table [TJ and compared it 
to the results obtained with the FE solver. A contour plot of vorticity, comparing the hybrid results to the 
FE ones is shown in Figure [I3| In Figure [l4j left, we plot the time evolution of the maximum of vorticity, 
showing the transfer of information between the two sub-domains of the hybrid solver. 

In Figure |14[ right, we compare the evolution of the vorticity maximum for different particle sizes, 
h = 5 x 10 -3 and ft = lx 10 -2 . As can be seen, for larger particle sizes it is possible to notice a strong 
decrease in the vorticity maximum as the dipole enters the Eulerian sub-domain. This produces a small 
reduction in the mean propagation speed of the dipole, as can be seen in Figure where the hybrid 
solution is slightly lagging behind the FE one. 


Table 1: Summary of the parameters used in the hybrid simulation of the Clercx-Bruneau dipole convection problem. 


Parameters 

Value 

Unit 

Description 

V 

1.6 x 1CT 3 

kgs -1 m _i 

Kinematic viscosity 

A 

1 

- 

Overlap ratio 

h 

5 x 1CT 3 

m 

Nominal blob spacing 

hgrid 

7 x 1CT 3 

m 

FE cell diameter 

N cells 

4 x 10 4 

- 

Number of mesh cells 

dbdry 

2h 

m 

Interpolation domain, 12/, offset from <912 e 

At L 

2.5 x 1CT 4 

s 

Lagrangian time step size 

AtE 

2.5 x 1CT 5 

s 

Eulerian time step size 
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Figure 13: Plot of the Clercx-Bruneau dipole at t = [0, 0.2, 0.4, 0.7]. The figure compares the hybrid simulation (top halves) 
against the FE only simulation (bottom halves). 




Figure 14: Evolution of the maximum of vorticity, max{cj}, from t = 0 to t = 0.7. The solutions are compared to FE [—, 
solid black] and VPM [- -, dashed black] benchmark simulations with characteristics identical to the Eulerian and Lagrangian 
components of the hybrid simulation. The figure depicts (left) the maximum vorticity in the Eulerian and Lagrangian sub- 
domains of the hybrid method for a blob spacing h = 0.005, and (right) the maximum vorticity of the hybrid method with 
nominal blob spacing h = 0.01 and h = 0.005. 


3.2. Dipole in bounded flow 

To further investigate the properties of the hybrid flow solver, namely its ability to capture the generation 
of vorticity at a solid wall, we applied it to the collision of the Clercx-Bruneau dipole with a solid wall, m 
A FE solution was validated against the results of Clercx and Bruneau, m , and used as reference. 

The setup of the hybrid domain is as shown in Figure [l5j The Eulerian sub-domain, £7^;, covers the 
near-wall region and the Lagrangian sub-domain domain extends to the complete fluid domain, which is 
bounded by the no-slip wall H w . The parameters used in this simulation are shown in Table [2] 
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Figure 15: The domain decomposition for the Clercx-Bruneau dipole collision problem, with the positive pole at = (x i, y\) = 
(0.1,0) and negative pole at p- = (x 2 ,V 2 ) = (—0.1,0). (Not to scale) 


Table 2: Summary of the parameters used in the hybrid simulation of the Clercx-Bruneau dipole wall collision problem. 


Parameters 

Value 

Unit 

Description 

V 

1.6 x 1CT 3 

kgs -1 m _i 

Kinematic viscosity 

A 

1 

- 

Overlap ratio 

h 

3 x 10" 3 

m 

Nominal blob spacing 

hgrid 

5 x 10" 3 to 1 x 10 -2 

m 

FE cell diameter 

N cells 

58 272 

- 

Number of mesh cells 

dbdry 

2 h 

m 

Interpolation domain, 0/ offset from 

dsurf 

3 h 

m 

Interpolation domain, 0/ offset from T, w 

A t L 

2.5 x 10” 4 

s 

Lagrangian time step size 

A t.E 

2.5 x 1CT 5 

s 

Eulerian time step size 


In Figure [16| we show contour plots of vorticity, comparing the hybrid results to the FE ones. We can see 
a good agreement between the two results. In Figure [17] we compare the evolution of vorticity maximum, 
energy, E, enstrophy, Q, and palinstrophy, P. As can be seen, the hybrid solver is capable of reproducing 
the results of the FE solver, with only small variations on the energy and palinstrophy. 
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Figure 16: Plot of the dipole at t = [0, 0.2, 0.4, 0.6, 0.8,1] (from left to right and top to bottom), comparing the hybrid simulation 
(left half) and FE only simulation (right half). 
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Figure 17: Variation in the fluid parameters from t = 0 to t = 1. The figure compares the hybrid results [—, solid blue] with 
the FE only [—, solid black] results. 
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3.3. Flow around cylinder, Re = 550 


An important aspect of flow simulations is the accurate calculation of forces, specifically lift and drag. 
Therefore, in this section we apply the hybrid solver to the flow around an impulsively started cylinder 
at Re = 550 and determine the forces acting on the cylinder. This test case problem has been extensively 
analysed in the literature, for example Koumoutaskos and Leonard, [35], and Rosenfeld et al., |S§], and these 
results will serve as a reference for the assessment of the hybrid solver, since one is a pure vortex particle 
solver and the other one is a pure Eulerian grid solver. 

The configuration of the hybrid domain is presented in Figure 18 As before, the Lagrangian domain, 
f II, covers the whole fluid domain and the Eulerian domain, C2#, is confined to a small region around the 
cylinder. The parameters used in this simulation are presented in Table [3] 



Figure 18: The domain decomposition for the impulsively started cylinder. The parameters of the domain are shown in Table 
[3] (Not to scale) 


Table 3: Summary of the parameters used in the hybrid simulation of the impulsively started cylinder at Re = 550. 


Parameters 

Value 

Unit 

Description 

C /00 

1.0e x 

ms -1 

Free stream velocity 

V 

3.6 x 10~ 3 

kgs -1 m _1 

Kinematic viscosity 

A 

1 

- 

Overlap ratio 

h 

8 x 10 -3 

m 

Nominal blob spacing 

hgrid 

8 x 10 -3 to 4 x 10 -2 

m 

FE cell diameter 

N cells 

32138 

- 

Number of mesh cells 

R 

1.0 

m 

Radius of the cylinder 

Rext 

1.5 

m 

Radius of the external boundary E^ 

dbdry 

0.21? 

m 

Interpolation domain, 0/ offset from E^ 

dsurf 

3 h 

m 

Interpolation domain, Qj offset from E„, 

A t.L 

3.0 x 10~ 3 

s 

Lagrangian time step size 

A tE 

1.0 x 10” 3 

s 

Eulerian time step size 


The contour plots of vorticity, comparing the hybrid results to the FE ones, are presented in Figure 19 


We can see that the two solvers give very similar results. Regarding the drag and lift, we can say that the 
hybrid solver is capable of reproducing both the FE results and the results of Koumoutsakos and Leonard, 
155]. see Figure 20 left. The differences exist mainly in the first instants of the simulation. It is important 
to highlight that, increasing h^dry from 0.11? to 0.21? improved considerably the results, see Figure [20] right. 
We have noted that this parameter is important and further research should be done in order to find its 
optimal value. A longer time simulation, t £ [0,40], was performed and the lift and drag compared to the 
results of Rosenfeld et al. [59], see Figure [ 2 T] The hybrid solver follows very well both the reference results 
of Rosenfeld et al. and of the FE simulation up to the onset of the vortex shedding, t ss 5s. After that, 
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all results stop having a good match but all remain within the same bounds and show similar frequency, as 
expected. 
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Figure 19: Plot of the vorticity field at t = [10,20,30,40], comparing the hybrid simulation (left) with the FE simulation 
(right). 




Figure 20: Left: Evolution of the drag coefficient during the initial stages t = 0 to t = 4 with total drag coefficient Cd (solid), 
pressure drag coefficient Cd pres (dashed) and friction drag coefficient Cdf ric (dotted). The figure compares results of hybrid 
simulation (blue), FE simulation (red) and reference data (black) of Koumoutsakos and Leonard m- Right: Comparison of 
total drag coefficient Cd for dydry £ {0.1 /i. 0.2 R \. 


Another aspect of the hybrid method, which is inherited from the vortex particle method, is its automatic 
adaptivity, where computational elements exist only in the support of vorticity, as opposed to standard grid 


solvers where the computational elements exist in the whole computational domain, see Figure 22 
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Figure 21: Evolution of the lift and drag coefficients from t = 0 to t = 40 with artificial perturbation m- The figure compares 
hybrid (blue), FE only (red), and the reference data (black) from Rosenfeld et al. Ell- 



Figure 22: Spatial distribution of computational elements for the hybrid method (left) and the FE method (right). 



3-4- Flow around stalled ellipsoid, Re = 5000 

The final test case consists in the simulation of the flow around a stalled ellipsoid at a Reynolds number 
Re = 5000. With this test case we aimed to assess the hybrid method’s capability to simulate flows with 
higher Reynolds number and to evaluate its performance in a situation where vortices exit and re-enter the 
Eulerian sub-domain. 

The configuration of the hybrid domain is presented in Figure |23[ Once more, the Lagrangian sub- 
domain, fi/,, covers the whole fluid domain and the Eulerian domain, fig, is restricted to the vicinity of the 
solid boundary. The parameters used in this simulation are presented in Table [4] 

Uco 



Figure 23: The domain decomposition for the stalled elliptical airfoil test case. The parameters of the domain are tabulated in 
Table [ 4 ] (Not to scale) 

The contour plots of vorticity, comparing the hybrid results to the FE ones, are presented in Figure [24} 
We can see that the two solvers give very similar results up to t = 3s and after that they start to diverge. 
This behaviour is expected due to the non-linear nature of the problem. Nevertheless we can see a compa¬ 
rable vortical structure. Regarding the drag and lift, we can observe that the hybrid solver is capable of 
reproducing very well the FE results up to t = 2s, see Figure [25] left. After that instant, which corresponds 
to the start of the vortex shedding, the results stop having a good match, but remain within the same 
bounds and show similar behaviour. 
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Table 4: Summary of the parameters used in the hybrid simulation of the stalled ellipsoid at Re = 5000. 


Parameters 

Value 

Unit 

Description 

E/oc 

1.0e x 

ms” 1 

Free stream velocity 

V 

2 x 10^ 4 

kg s” 1 m” 1 

Kinematic viscosity 

A 

1 

- 

Overlap ratio 

h 

1.67 x 10" 3 

m 

Nominal blob spacing 

fogrid 

1.4 x 1CT 3 to 8 x 10’ 3 

m 

FE cell diameter 

N cells 

118196 

- 

Number of mesh cells 

L 

1.0 

m 

Chord length 

Lext 

1.0 

m 

Length of the external boundary 

w 

0.12 

m 

Maximum thickness 

w ext 

1.0 

m 

Maximum thickness of the external boundary 

Rext 

1.5 

m 

Radius of the external boundary E^ 

dbdry 

0.1L 

m 

Interpolation domain, flj offset from E^ 

dsurf 

3 h 

m 

Interpolation domain, flj offset from E u , 

A t L 

3.0 x 10" 3 

s 

Lagrangian time step size 

A tE 

1.0 x 10" 3 

s 

Eulerian time step size 
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Figure 24: Plot of the vorticity field for t E {1-0, 2.0, 3.0, 3.5, 4.0}, comparing the hybrid simulation (left) with the FE simulation 
(right). 


23 


































































































































































Figure 25: Evolution of the lift and drag coefficient from t = 0 up to t = 10. The figure compares the hybrid results with FE 
simulation. 
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4. Conclusions and outlook 


In this work we have presented an efficient hybrid Eulerian-Lagrangian flow solver that does not rely 
on the Schwarz iteration method and is capable of exactly preserving the total circulation, improving the 
results of similar approaches. 

By dividing the computational domain into two types of sub-domains, Lagrangian based and Eulerian 
based, it is possible to use the most suitable solver in each region of the flow domain. This formulation 
allows for the use of small and highly resolved Eulerian solvers in the near wall regions. These solvers can 
efficiently capture the generation of vorticity at the solid boundaries, potentially using wall functions and 
other advanced methods. In the wake region, the Lagrangian solver accurately models the development of 
the wake due to its practically absent numerical diffusion. 

This approach paves the way to the development of complex simulations of wind turbines since each 
solid object can be individually modelled by a disjoint Eulerian subdomain embedded in the Lagrangian 
one. This greatly simplifies the parallelizability of large flow simulations. 

The future developments of this work are the extension to moving and deformable objects, extension 
to large Reynolds number flows (including the coupling of turbulence between the two sub-domains) and 
extension to three dimensional flows. 


Acknowledgements 

This work is made possible in part due to funds from Technology Foundation STW - Veni grant, and 
Future Emerging Technologies FP7-Deepwind project and the U.S. Department of Energy. 


References 

[1] L. A. Barba, Vortex method for computing high-Reynolds number flows: Increased accuracy with a fully mesh-less 
formulation, Ph.D. thesis, California Institute of Technology (2004). 

[2] L. A. Barba, A. Leonard, C. B. Allen, Advances in viscous vortex methods—meshless spatial adaption based on radial 
basis function interpolation, International Journal for Numerical Methods in Fluids 47 (5) (2005) 387—421. 

[3] L. A. Barba, R. Yokota, How Will the Fast Multipole Method Fare in the Exascale Era?, SIAM News 46 (6). 

[4] J. Barnes, P. Hut, A hierarchical O(nlogn) force-calculation algorithm, Nature 324 (6096) (1986) 446-449. 

[5] J. Beale, On the Accuracy of Vortex Methods at Large Times, in: Computational Fluid Dynamics and Reacting Gas 
Flows SE - 2, vol. 12, Springer New York, 1988, pp. 19-32. 

[61 J. Beale, A. Maida, High order accurate vortex methods with explicit velocity kernels, Journal of Computational Physics 
58 (2) (1985) 188-208. 

[71 J. Beale, A. Maida, High order accurate vortex methods with explicit velocity kernels, Journal of Computational Physics 
58 (2) (1985) 188-208. 

[8] J. T. Beale, A. Majda, Vortex methods. II. Higher order accuracy in two and three dimensions, Mathematics of Compu¬ 
tation 39 (159) (1982) 29-52. 

[9] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, vol. 15, Springer-Verlag, New York, 1991. 

[10] R. Chatelin, P. Poncet, Hybrid grid-particle methods and Penalization: A Sherman-Morrison-Woodbury approach to 
compute 3D viscous flows using FFT, Journal of Computational Physics 269 (2014) 314—328. 

[11] A. J. Chorin, Numerical solution of the Navier-Stokes equations, Mathematics of Computation 22 (104) (1968) 745—762. 

[12] A. J. Chorin, Numerical study of slightly viscous flow, Journal of Fluid Mechanics 57 (04) (1973) 785-796. 

[13] H. Clercx, C.-H. Bruneau, The normal and oblique collision of a dipole with a no-slip boundary, Computers & Fluids 
35 (3) (2006) 245-279. 

[14] G. H. Cottet, Methodes particulaires pour l’equation d’Euler dans le plan, Ph.D. thesis, Universite Paris VI, Paris (1982). 

[15] G. H. Cottet, Particle-grid domain decomposition methods for the Navier-Stokes equations in exterior domains, in: Vortex 

Dynamics and Vortex Methods, American Mathematical Society, 1991, pp. 103—117. 

[16] G.-H. Cottet, A vorticity creation algorithm for the Navier-Stokes equations in arbitrary domain, in: Navier-Stokes 
equations and related non-linear problems, 1994. 

[17] G.-H. Cottet, P. D. Koumoutsakos, Vortex Methods: Theory and Practice, Cambridge University Press, 2000. 

[18] G.-H. Cottet, M.-L. Ould Salihi, M. El Hamroui, Multi-purpose regridding in vortex methods, ESAIM: Proceedings 7 

(1999) 94-103. 

[19] G. Daeninck, Developments in hybrid approaches : Vortex method with known separation location, Ph.D. thesis, Universite 
Catolique de Louvain (2006). 

[20] P. Degond, S. Mas-Gallic, The Weighted Particle Method for Convection-Diffusion Equations Part 1: The Case of an 
Isotropic Viscosity, Mathematics of Computation 53 (1989) 485—507. 


25 



[21] S. Dong, G. E. Karniadakis, C. Chryssostomidis, A robust and accurate outflow boundary condition for incompressible 
flow simulations on severely-truncated unbounded domains, Journal of Computational Physics 261 (2014) 83-105. 

[22] S. Engblom, On well-separated sets and fast multipole methods, Applied Numerical Mathematics 61 (10) (2011) 1096-1102. 

[23] K. Goda, A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows, 
Journal of Computational Physics 30 (1) (1979) 76-95. 

[24] A. Golas, R. Narain, J. Sewall, P. Krajcevski, P. Dubey, M. Lin, Large-scale fluid simulation using velocity-vorticity 
domain decomposition, ACM Transactions on Graphics 31 (6) (2012) 1. 

[25] A. Goude, S. Engblom, Adaptive fast multipole methods on the GPU, The Journal of Supercomputing 63 (3) (2012) 
897-918. 

[26] L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, Journal of Computational Physics 73 (2) (1987) 
325-348. 

[27] P. M. Gresho, Incompressible fluid dynamics: Some fundamental formulation issues, Annual Review of Fluid Mechanics 
23 (1) (1991) 413-453. 

[28] J. Guermond, P. Minev, J. Shen, An overview of projection methods for incompressible flows, Computer Methods in 
Applied Mechanics and Engineering 195 (44-47) (2006) 6011-6045. 

[29] J.-L. Guermond, S. Huberson, W.-Z. Shen, Simulation of 2D External Viscous Flows by Means of a Domain Decomposition 
Method, Journal of Computational Physics 108 (2) (1993) 343-352. 

[30] J.-L. Guermond, W.-Z. Shen, A domain decomposition method for simulating 2D external viscous flows, in: Domain 
Decomposition Methods in Science and Engineering: The Sixth International Conference on Domain Decomposition, 
American Mathematical Soc., 1994, pp. 463-467. 

[31] O. H. Hald, Convergence of Vortex Methods for Eulers Equations. II, SIAM Journal on Numerical Analysis 16 (5) (1979) 
726-755. 

[32] J. G. Heywood, R. Rannacher, S. Turek, Artificial boundaries and flux and pressure conditions for the incompressible 
Navier-Stokes equations, International Journal for Numerical Methods in Fluids 22 (5) (1996) 325—352. 

[33] L. H. Howell, J. B. Bell, An adaptive mesh projection method for viscous incompressible flow, SIAM Journal on Scientific 
Computing 18 (4) (1997) 996-1013. 

[34] S. G. Huberson, S. G. Voutsinas, Particles and grid, Computers Sz Fluids 31 (4-7) (2002) 607-625. 

[35] P. Koumoutsakos, Direct numerical simulations of unsteady separated flows using vortex methods, Ph.D. thesis, California 
Institute of Technology (1993). 

[36] P. Koumoutsakos, Inviscid Axisymmetrization of an Elliptical Vortex, Journal of Computational Physics 138 (2) (1997) 
821-857. 

[37] P. Koumoutsakos, A. Leonard, Improved boundary integral method for inviscid boundary condition applications, AIAA 
Journal 31 (2) (1993) 401-404. 

[38] P. Koumoutsakos, A. Leonard, High-resolution simulations of the flow around an impulsively started cylinder using vortex 
methods, Journal of Fluid Mechanics 296 (1) (1995) 1-38. 

[39] P. Koumoutsakos, A. Leonard, F. Pepin, Boundary Conditions for Viscous Vortex Methods, Journal of Computational 
Physics 113 (1) (1994) 52-61. 

[40] P. Koumoutsakos, D. Shiels, Simulations of the viscous flow normal to an impulsively started and uniformly accelerated 
flat plate, Journal of Fluid Mechanics 328 (1996) 177-227. 

[41] I. Lashuk, G. Biros, A. Chandramowlishwaran, H. Langston, T.-A. Nguyen, R. Sampath, A. Shringarpure, R. Vuduc, 
L. Ying, D. Zorin, A massively parallel adaptive fast multipole method on heterogeneous architectures, Communications 
of the ACM 55 (5) (2012) 101. 

[42] Y. Lecointe, J. Piquet, On the use of several compact methods for the study of unsteady incompressible viscous flow round 
a circular cylinder, Computers Sz Fluids 12 (4) (1984) 255-280. 

[43] A. Leonard, Computing Three-Dimensional Incompressible Flows with Vortex Elements, Annual Review of Fluid Mechan¬ 
ics 17 (1985) 523-559. 

[44] M. Lesoinne, C. Farhat, Geometric conservation laws for flow problems with moving boundaries and deformable meshes, 
and their impact on aeroelastic computations, Computer Methods in Applied Mechanics and Engineering 134 (1-2) (1996) 
71-90. 

[45] M. J. Lighthill, Introduction: boundary layer theory, in: Laminar boundary layers, Oxford University Press, 1963, pp. 
46-113. 

[46] A. Logg, G. N. Wells, DOLFIN, ACM Transactions on Mathematical Software 37 (2) (2010) 1-28. 

[47] L. Manickathan, Hybrid Eulerian-Lagrangian vortex particle method, Master’s thesis, Delft University of Technology 
(2015). 

[48] J. Monaghan, Extrapolating B-splines for interpolation, Journal of Computational Physics 60 (2) (1985) 253-262. 

[49] G. Morgenthal, J. Walther, An immersed interface method for the Vortex-In-Cell algorithm, Computers Sz Structures 
85 (11-14) (2007) 712-726. 

[50] H. O. Nordmark, Rezoning for higher order vortex methods, Journal of Computational Physics 97 (2) (1991) 366-397. 

[51] M. L. Ould-Salihi, G. H. Cottet, M. El Hamraoui, Blending Finite-Difference and Vortex Methods for Incompressible Flow 
Computations, SIAM Journal on Scientific Computing 22 (5) (2001) 1655—1674. 

[52] G. S. Oxley, A 2-D Hybrid Euler-Compressible Vortex Particle Method For Transonic Rotorcraft Flows, Ph.D. thesis, 
Carleton University (2009). 

[53] G. Papadakis, S. G. Voutsinas, In view of accelerating CFD simulations through coupling with vortex particle approxi¬ 
mations, Journal of Physics: Conference Series 524 (1) (2014) 012126. 

[54] P. Ploumhans, G. Winckelmans, Vortex Methods for High-Resolution Simulations of Viscous Flow Past Bluff Bodies of 


26 



General Geometry, Journal of Computational Physics 165 (2) (2000) 354-406. 

[55] P. Ploumhans, G. Winckelmans, J. Salmon, A. Leonard, M. Warren, Vortex Methods for Direct Numerical Simulation of 
Three-Dimensional Bluff Body Flows: Application to the Sphere at Re=300, 500, and 1000, Journal of Computational 
Physics 178 (2) (2002) 427-463. 

[56] J. T. J. Rasmussen, G.-H. Cottet, J. H. Walther, A multiresolution remeshed Vortex-In-Cell algorithm using patches, 
Journal of Computational Physics 230 (17) (2011) 6742-6755. 

[57] P. Raviart, An analysis of particle methods, in: Numerical Methods in Fluid Dynamics, vol. 1127, Springer Berlin 
Heidelberg, 1985, pp. 243-324. 

[58] F. Renac, S. Gerald, C. Marmignon, F. Coquel, Fast time implicit-explicit discontinuous Galerkin method for the com¬ 
pressible Navier-Stokes equations, Journal of Computational Physics 251 (2013) 272-291. 

[59] M. Rosenfeld, D. Kwak, M. Vinokur, A fractional step solution method for the unsteady incompressible Navier-Stokes 
equations in generalized coordinate systems, Journal of Computational Physics 94 (1) (1991) 102-137. 

[60] L. Rosenhead, The Spread of Vorticity in the Wake Behind a Cylinder, Proceedings of the Royal Society of London A: 
Mathematical, Physical and Engineering Sciences 127 (806) (1930) 590-612. 

[61] R. L. Sani, P. M. Gresho, Resume and remarks on the open boundary condition minisymposium, International Journal 
for Numerical Methods in Fluids 18 (10) (1994) 983-1008. 

[62] I. F. Sbalzarini, J. H. Walther, M. Bergdorf, S. E. Hieber, E. M. Kotsalis, P. Koumoutsakos, PPM - A highly efficient 
parallel particle-mesh library for the simulation of continuum systems, Journal of Computational Physics 215 (2) (2006) 
566-588. 

[63] S. Shankar, L. van Dommelen, A New Diffusion Procedure for Vortex Methods, Journal of Computational Physics 127 (1) 
(1996) 88-109. 

[64] D. Shiels, Simulation of controlled bluff body flow with a viscous vortex method, Ph.D. thesis, California Institute of 
Technology (1998). 

[65] D. R. Speck, Generalized Algebraic Kernels and Multipole Expansions for massively parallel Vortex Particle Methods, 
Ph.D. thesis, Bergische Universitat Wuppertal (2011). 

[66] M. Stock, A. Gharakhani, C. Stone, Modeling Rotor Wakes with a Hybrid OVERFLOW-Vortex Method on a GPU Cluster, 
in: 28th AIAA Applied Aerodynamics Conference, 2010. 

[67] C. Taylor, P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Computers & 
Fluids 1 (1) (1973) 73-100. 

[68] T. E. Tezduyar, Finite Element Methods for Flow Problems with Moving Boundaries and Interfaces, Archives of Compu¬ 
tational Methods in Engineering 8 (2) (2001) 83-130. 

[69] D. Wee, A. F. Ghoniem, Modified interpolation kernels for treating diffusion and remeshing in vortex methods, Journal 
of Computational Physics 213 (1) (2006) 239-263. 

[70] G. Winckelmans, A. Leonard, Contributions to Vortex Particle Methods for the Computation of Three-Dimensional 
Incompressible Unsteady Flows, Journal of Computational Physics 109 (2) (1993) 247-273. 

[71] G. S. Winckelmans, Topics in vortex methods for the computation of three- and two-dimensional incompressible unsteady 
flows, Ph.D. thesis, California Institute of Technology, USA (1989). 

[72] G. S. Winckelmans, Vortex methods, in: Encyclopedia of Computational Mechanics, vol. 3, John Wiley Sz Sons, 2004, pp. 
129-153. 

[73] R. Yokota, L. A. Barba, A tuned and scalable fast multipole method as a preeminent algorithm for exascale systems, 
International Journal of High Performance Computing Applications 26 (4) (2012) 337-346. 


27 



